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Abstract 

The vertical distribution of phytoplankton is of fundamental importance for the dynamics and structure of aquatic 
communities. Here, using an advection-reaction-diffusion model, we investigate the distribution and competition of 
phytoplankton species in a water column, in which inverse resource gradients of light and a nutrient can limit growth 
of the biomass. This problem poses a challenge for ecologists, as the location of a production layer is not fixed, but 
rather depends on many internal parameters and environmental factors. In particular, we study the influence of an 
upper mixed layer (UML) in this system and show that it leads to a variety of dynamic efl'ects: (i) Our model predicts 
alternative density profiles with a maximum of biomass either within or below the UML, thereby the system may be 
bistable or the relaxation from an unstable state may require a long-lasting transition, (ii) Reduced mixing in the deep 
layer can induce oscillations of the biomass; we show that a UML can sustain these oscillations even if the diff'usivity 
is less than the critical mixing for a sinking phytoplankton population, (iii) A UML can strongly modify the outcome 
of competition between diff'erent phytoplankton species, yielding bistability both in the spatial distribution and in the 
species composition, (iv) A light limited species can obtain a competitive advantage if the diff'usivity in the deep 
layers is reduced below a critical value. This yields a subtle competitive exclusion eff'ect, where the oscillatory states 
in the deep layers are displaced by steady solutions in the UML. Finally, we present a novel graphical approach for 
deducing the competition outcome and for the analysis of the role of a UML in aquatic systems. 
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1. Introduction 

The survival and competition of species in an hetero- 
geneous environment has fascinated ecologists for long 
times (see e.g. Holmes et al. 1994; Tilman and Kareiva 
1997; Huisman et al. 1999; Neuhauser 2001). In many 
systems the spatial diversity of natural populations orig- 
inates mainly from some underlying abiotic heterogene- 
ity of the environment. If growth conditions vary be- 
tween diff'erent locations this spatial variation should 
be reflected in the density distribution of natural pop- 
ulations. After the seminal papers by Skellam (1951) 
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and by Kierstead and Slobodkin (1953), the dynam- 
ics of such systems have often been analyzed in terms 
of favorable and unfavorable patches (see e.g. Okubo 
and Levin 2001; CantreU and Cosner 2001; Birch et al. 
2007). This approach assumes that space can be divided 
into regions of positive and negative net growth, be- 
tween which the organisms are transported by difl'usion 
and advection. These suggestions, although realistic in 
many situations, do not hold in some resource-consumer 
systems, in which the size, the form and even the loca- 
tion of the species' favorable patch may vary, reflecting 
external and internal perturbations. 

An important example is provided by the dynamics 
of phytoplankton in an incompletely mixed water col- 
umn (see e.g. Jamart et al. 1977; Klausmeier and Litch- 
man 2001; Beckmann and Hense 2007). In many re- 
gions of the world's ocean the form of vertical phy- 
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toplankton profiles is determined mainly by two fac- 
tors. First, the reduction of the light intensity with depth 
makes deep layers unfavorable for photosynthetic phy- 
toplankton species. Second, an opposing gradient of nu- 
trients can often maintain a positive net production only 
in deep subsurface layers. These deep states, known 
as deep chlorophyll maxima (DCM), constitute one of 
the most striking characteristics of nutrient poor wa- 
ters in ocean ecosystems and freshwater lakes (Steele 
and Yentsch 1960; Anderson 1969; Abbott et al. 1984; 
Cullen 1982; Tittel et al. 2003). The production layer in 
such systems is highly variable and dynamic. First it is 
species specific: less nutrient limited species will have 
a positive production rate close to the surface, whereas 
smaller light requirements will result in subsurface pro- 
ductivity (see e.g. Venrick 1993). Second, as each 
species shades light and consumes nutrients it acts as 
an ecosystem-engineer, modifying its abiotic environ- 
ment (Jones et al. 1994). As a consequence, the posi- 
tion of the favorable layer for a phytoplankton species 
depends on the current abundance and distribution of its 
own biomass, as well as that of all other species. 

Theoretical models have been a useful approach to 
describe nutrient limited phytoplankton growth in con- 
stant and seasonally driven environments (Huppert et al. 
2002, 2005). The dynamics, competition and structur- 
ing within a vertical water column have been investi- 
gated in a series of modeling investigations (Radach and 
Maier-Reimer 1975; Shigesada and Okubo 1981; Varela 
et al. 1992; Huisman and Weissing 1995; Huisman et 
al. 1999; Klausmeier and Litchman 2001; Diehl 2002; 
Hodges and Rudnick 2004; Huisman et al. 2006; Beck- 
mann and Hense 2007). It was shown that a given set of 
parameters may lead to either a surface or a deep chloro- 
phyll maximum, whereby the location of the maximum 
is entirely determined by the environmental conditions. 
These model results are in agreement with many field 
studies (e.g., Aristegui et al. 2003; Matondkar et al. 
2005; Weston et al. 2005). However, in a few recent 
investigations either surface or deep maxima were ob- 
served under almost the same conditions (Venrick 1993; 
Holm-Hansen and Hewes 2004). These observations 
suggest the existence of bistability in the spatial phy- 
toplankton configuration, however they cannot easily 
be reproduced in current models and may indicate that 
some important mechanisms, which are contributing to 
the spatial organization of aquatic communities, are still 
not well understood. 

One major ingredient determining whether a species 
can establish close to the surface is the presence of an 
upper mixed layer (UML). UMLs commonly occur in 
oceans and lakes due to mechanical perturbation of the 



surface waters (e.g., by wind, waves, storms) and are 
characterized by strong turbulent mixing up to a depth 
of about 30 m to 200 m or more (Denser 1987; Venrick 
1993; Law et al. 2000). To our knowledge the first study 
to demonstrate the strong influence of a UML on the 
spatial configuration of phytoplankton was performed 
by Yoshiyama and Nakajima (2002, 2006), who consid- 
ered a single species model where the water column was 
divided into an infinitely mixed UML and poorly mixed 
lower layers, with a small diff'usivity across the sepa- 
rating boundary layer (thermocline). The model out- 
come were vertical phytoplankton patterns with a sharp 
edge at the thermocline, with the possibility of bistabil- 
ity in the phytoplankton distribution, characterized by 
phase transitions between biomass maxima in the deep 
and upper layers. However, these articles did not ana- 
lyze such important questions as the influence of finite 
mixing in the upper layer, the competition of species 
occupying diff'erent layers, or the dynamics of the sys- 
tem when deep phytoplankton maxima are not station- 
ary, and thus many important questions about the role of 
the UML for structuring phytoplankton configurations 
still remain open. 

In this study, we consider a mathematical model for 
the growth of a nutrient and light limited phytoplank- 
ton community in a vertical incompletely mixed water 
column. A UML is introduced by smoothly varying the 
strength of diffusivity with vertical depth, so that dis- 
continuities are avoided. We find that a UML amelio- 
rates the growth conditions close to the surface and can 
have drastic eff'ects because organisms close to the sur- 
face occlude light and prevent growth in all deeper lay- 
ers. Thus the presence or absence of a UML turns out 
to be a major factor controlling the vertical distribution 
and competition outcome in the whole water column. 

In particular, the presence of a UML yields the fol- 
lowing key model outcomes: First, the spatial density 
profile can become bistable with vertical maxima either 
close to the surface or in deep layers. Close to the bista- 
bility range dynamics are characterized by time scale 
separation and the relaxation from an unstable state may 
require a long-lasting transition. Second, in the pres- 
ence of a UML the persistence threshold of sinking phy- 
toplankton disappears. Thus, our model provides a new 
mechanism for the persistence of a population in a flow 
(Hershey et al. 1993; Huisman et al. 2002; Ryabov 
and Blasius 2008). Third, a UML can strongly modify 
the competition outcome between diflferent phytoplank- 
ton species by providing a vertical niche for species 
which are better adapted to the conditions close to the 
surface. As a result, in the two species model, we ob- 
serve bistability both in the spatial distribution and in 
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the species composition. Fourth, we show that a light 
Hmited species can obtain a competitive advantage if the 
difFusivity in the deep layers is reduced below a critical 
value. In this case we identify a dynamic competitive 
exclusion effect, where the species composition is deter- 
mined by the dynamic state such that oscillatory states 
in the deep layers are displaced by steady solutions in 
the UML. Finally, we present a graphical approach for 
analysing the mechanism of resource competition in a 
spatially extended system. 

2. Model 

Our model describes the dynamics of n phytoplank- 
ton species, competing for resources in a vertical wa- 
ter column of depth Zg (in the following we consider 
only the cases n = I and n = 2). We study a bottom- 
up control of phytoplankton, which means that we did 
not consider the influence of zooplankton and higher 
trophic levels. Let Pi(z, t) denote the density of species 
/ at time t and depth z. Assume that there are two limit- 
ing factors: the concentration of a nutrient, N{z, t), and 
the light intensity I(z, t), both of which are a function of 
the vertical position z. Coupling of the nutrient and the 
phytoplankton dynamics leads to the following system 
of reaction-advection-difl'usion equations (Radach and 
Maier-Reimer 1975; Jamart et al. 1977; Klausmeier and 
Litchman 2001; Huisman et al. 2006) 

dPi 



where yUmax is the maximum growth rate, and H^^ and 
H^^^ are the corresponding half saturation constants. 
Varying and /// allows to model, for instance, a 
species which is better adapted for light (a smaller Hj) 
or for the nutrient (a smaller Hn). 

Light dissipates with depth as it is absorbed by the 
biomass, water, clay particles and many other absorbing 
substances. Assume that the light intensity decreases 
exponentially according to Lamber-Beer's law (see e.g. 
Shigesada and Okubo 1981, Kirk 1994) 
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where jdi(N, I) describes the local growth rate of species 
/, m is the mortality, v is the phytoplankton sinking ve- 
locity, D(z) is the depth dependent turbulent difl'usivity, 
a is the nutrient content of a phytoplankton cell and s is 
the phytoplankton recycling coefficient. 

Equations ([B and Q are coupled by means of the 
growth rate jdi(N, I) which depends on the local resource 
availability and also controls the nutrient uptake. As- 
suming that the limitation of growth follows the Monod 
kinetics (e.g. Turpin 1988) and both resources are es- 
sential (von Liebig's law of minimum), we obtain 
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I(z) = lin exp 
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where //„ is the incident light intensity, K^g is the tur- 
bidity of the water without biomass, and k is the phyto- 
plankton light-absorption coefficient. 

To describe the water column stratification, we as- 
sume that the diff'usivity D(z) takes a large value Du in 
the upper layer and a much smaller value Do in the deep 
layers. A gradual transition from one area to another can 
be written in terms of a generalized Fermi function 



D(z) = Dd + 



Du-Dd 



(5) 



where Zmix describes the depth of the UML and the pa- 
rameter w characterizes the width of the transient layer. 
In all numerical experiments we chose Du = 50 cm^/s, 
modeling well-mixed waters in the UML, and Do = 
0.1. ..1.0 cm^/s for the lower layers (Lewis et al. 1986; 
Saggio and Imberger 2001 ; Smyth et al. 2001 ; Finnigan 
et al. 2002). The model parameters were chosen to de- 
scribe clear ocean water (Huisman et al. 2006) and can 
be found in Table [T] 

We simulate the model up to a depth of z = Zb and 
assume impenetrable boundaries at the surface and at 
the bottom for the phytoplankton 
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Since at the bottom of a sufficiently deep water column 
(we assume Zb = 300 m) the phytoplankton biomass 
vanishes together with flux, one can also use the Dirich- 
let boundary condition at the bottom P/(Zfi) = and this 
will not influence our results. For the nutrient distribu- 
tion we assume an impenetrable surface and a constant 
concentration at the bottom 



= 0, N(Zb) = Nb . 



(7) 



z=0 



As initial condition we assumed a small uniformly 
distributed concentration of phytoplankton (P(z,0) < 



3 THE SINGLE SPECIES MODEL 



A. Ryabov et al./ J. Theor. Biol. 263 (2010) 120-133 



1 cell m"^), whereas for the nutrient we used two differ- 
ent initial profiles, describing a nutrient saturated water 
column (N(z, 0) = Nb) and a nutrient depleted upper 
layer (N(z, 0) = if z < Z^ix)- We have also explored 
the influence of other initial conditions, however we did 
not find any new solutions beside the ones described. 
Checking the stability of solutions, in each case we sim- 
ulated dynamics for 50,000 system days (approximately 
130 years). 

The model was integrated using a backward diff'er- 
ence method, based on the finite volume scheme (Pham 
Thi et al. 2005). For the numerical solution we have 
discretized all variables on a grid which consisted of 
600 points. Diff'usion terms were approximated by 
a second order central discretization scheme, the ad- 
vection term was represented by a third-order upwind 
biased formula, integration was made via the trape- 
zoidal rule. The resulting system of ordinary diff'er- 
ential equations was solved by the CVODE package 
(littp://www.netlib.org/ode). For model validation we 
have compared our simulation results with already pub- 
lished results (Huisman et al. 2006) and further verified 
that the results remain unchanged if we double the num- 
ber of points in the grid. Furthermore, in some limiting 
cases it was possible to compare our simulation results 
with analytic solutions. 

3. The single species model 

We first concentrate on the dynamics of a single 
species population and describe the formation of a DCM 
in a water column without a UML. Suppose that we start 
with an initially nutrient rich system (N(z, t) = Nb). 
Thereby the nutrient limitation is negligible and we ob- 
serve a rapid formation of a transient phytoplankton 
maximum close to the surface (Fig. [T]\). This phyto- 
plankton profile P(z, t) is, however, not stable. With the 
depletion of the nutrient in the surface layer the pro- 
duction layer, i.e. the layer where jdiN, I) > m, shifts 
downwards (see also Fig.[2l\), until the system reaches 
a stable DCM configuration (Fig. [T^). In this equilib- 
rium the upward flux of nutrient compensates the nutri- 
ent consumption and any further sinking of the produc- 
tion layer is balanced by light limitation (Klausmeier 
and Litchman 2001). 

Note that the spatio-temporal evolution of the con- 
centration profile P(z, t) depends on the growth condi- 
tions at all vertical positions. Due to the water turbid- 
ity and phytoplankton self- shading the light intensity 
/(z, t) is reduced in deeper layers, thereby increasing 
the limitation by light (dashed lines in Fig.[T]). In con- 
trast, the nutrient concentration N{z, t) is close to zero 



within the bulk of phytoplankton biomass and above 
it, and increases almost linearly with depth below the 
phytoplankton peak (see dot-dashed lines in Fig. [T] for 
the nutrient limitation of growth). This shaping of the 
spatial dependence of the growth limiting factors self- 
consistently depends on the full phytoplankton density 
profile P{z, t), a fact which makes the problem very hard 
to understand without mathematical simulation. 

A rough insight into the time evolution of P(zj) 
can be gained by considering the centers of biomass 
Zm = f zPdzl f Pdz and of phytoplankton net produc- 
tion Zg = J zgPdzl f gPdz (black and gray arrows in 
Fig.[T]). Here g(z) describes the net phytoplankton pro- 
duction which includes growth, loss, sinking and mix- 
ing, i.e., the product gP equals the right-hand- side of 
Eq. (O. In an incompletely mixed water column with- 
out a UML the position of the center of mass, Z^, fol- 
lows that of Zg. Thereby the phytoplankton production 
around Zg shifts the mass center Z^, changing the local 
nutrient consumption and the light absorption, which in 
turn has an influence again on the position of the growth 
center Zg. In this feedback-loop, the system eventu- 
ally reaches a stable equilibrium configuration where 
both centers coincide, Zg = Zm, giving rise to a DCM 
(Fig.IB). 

Now suppose that there is strong mixing in the upper 
layer. If the bulk biomass is located sufficiently deep, 
then the mixing in the upper layer has practically no 
eff'ect and an identical DCM can persist (Fig. [T]C) in- 
dependent of whether or not a UML is present. Note 
that we always assume that the depth of the UML is 
smaller than the compensation depth (i.e., a depth at 
which yu(/) = min the absence of biomass, see e.g. Sver- 
drup 1953, Huisman et al. 1999b). 

By contrast, if the phytoplankton biomass is initially 
located close to the surface, it will be almost uniformly 
distributed within the UML. The position of is then 
fixed approximately in the middle of the UML and is 
almost uncoupled from Zg (Fig.[TJ)). Therefore a grad- 
ual shift of the center of mass into the deep layers is 
no longer possible and the transition to a DCM can 
only take place if the light intensity below the UML 
is sufficiently large to provide positive net growth in 
the deep layers - otherwise the phytoplankton remains 
trapped in the UML. We denote this stable configuration 
of a nearly uniform phytoplankton profile in the UML 
as an upper chlorophyll maximum (UCM, note that all 
acronyms are listed in Table O. Figs.[T]C and[Tt) show 
that in the presence of a UML, depending on the initial 
conditions, the system can undergo two very diff'erent 
spatial configurations of either a deep or an upper phy- 
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Figure 1 : Typical vertical phytoplankton profiles in the single species model for a system without a UML (top) and with a UML (bottom, depth 
Z,nix indicated by the horizontal dashed line). Plotted are the density of phytoplankton as function of depth, P(z) (black solid line), and the growth 
limiting terms with respect to light, 1/(1 + Hi) (gray dashed line), and to the nutrient, N/(N + H]\j) (gray dot-dashed line); positive growth is possible 
where both regulating terms are larger than m/yUmax (level of zero net growth indicated by the vertical dotted line). Black and gray arrows show 
the centers of biomass and net production Zg, respectively. (A) without a UML, a non-stable phytoplankton maximum close to the surface is 
driven downwards (indicated by the arrows) and evolves to the stable state (B). Under the same conditions in a system with a UML, we observe 
two alternative stable configurations: (C) a phytoplankton profile with a maximum in the deep layers (DCM), or (D) a profile with a maximum in 
the upper layer (UCM). 



toplankton maximum (Yoshiyama and Nakajima 2002). 

Fig. [2l depicts the typical spatio-temporal evolution of 
the phytoplankton density. Without a UML (Fig. [JK), 
an initial phytoplankton maximum at the surface slowly 
moves downward until the distribution converges to a 
stable DCM equilibrium. However, as is shown in 
Fig. [2^, DCM's are not necessarily steady states. For 
small values of diffusivity Do the sinking of biomass 
may destabilize the density profile, yielding sustained 
regular or chaotic oscillations of biomass in the deep 
layer (Huisman et al. 2006). Moreover, if the diff'usivity 
is lower than the minimal persistence threshold 

- 

4(/j(NBJin)-m) 

the sinking phytoplankton population cannot survive 
and goes extinct (Riley et al. 1949; Shigesada and 



Okubo 1981). Our numerical simulations show that in 
a system without a UML these three model outcomes 
(i.e., stationary DCM, oscillating DCM or extinction) 
are the only possible long-term solutions. Further, these 
solutions are globally attracting, which means that they 
are reached from any initial condition. 

In the presence of a UML these solutions can be 
found as well, however the dynamics may be more 
complicated. Most notably, as already mentioned, the 
system may be bistable: under the same parameters 
as in Fig. [2l\, an initially nutrient saturated water col- 
umn gives rise to a UCM (Fig. WP), whereas an ini- 
tially nutrient depleted water column leads to a DCM 
(Fig.[2t)). Note that we observe bistability only in a cer- 
tain parameter region in which a DCM is not aff'ected 
by the upper layer and a UCM contains enough biomass 
to limit growth in the deep layers. The diff'usivity in 
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Figure 2: Spatio-temporal evolution of typical phytoplankton profiles. Plotted is the phytoplankton density t) in color-coding (xlO^ cells m~^) 
for different values of Do and initial conditions, in a system without a UML (A-B) and in the presence of a UML (C-H). Black lines track the 
evolution of the center of biomass, Z^, white dashed lines indicate the depth of the UML. Dynamics without a UML: (A) gradual evolution of a 
DCM, Dd = 0.3 cm^/s, (B) oscillations of the biomass for small mixing, Do = 0.12 cm^/s. Bistability (with a UML): depending on the initial 
conditions either (C) a stable UCM or (D) a stable DCM is formed (value of as in A). Transition from an unstable state: for any initial condition 
(E) only a DCM is stable. Do = 0.2 cm^/s, or (F) only a UCM is stable. Do = 0.4 cm^/s, however the transient process may last a long time. 
Oscillations of the biomass: (G) oscillations are not affected by a UML, Do = 0.12 cm^/s (value of Do as in B), or (H) oscillations are induced by 
a UML, Dd = 0.04 cm^/s (value of Do < D^in = 0.0408 cm^/s). 
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the deep layers, Do, is a suitable bifurcation parame- 
ter as it controls the nutrient flux from the bottom and 
an increase of Do rises the level of the DCM (Klaus- 
meier and Litchman 2001). Thus decreasing Dd 
obtain steady (Figs. [2^) or oscillatory DCM's ([2lGr and 
[2^1), whereas larger values of Dd yield UCM solutions 
(Fig.Ef). 

Close to the bistability range we observe dynamics 
at very difl'erent time scales. Consider Fig. [2^ which 
shows the transition from an unstable upper maximum 
to a stable DCM. Here, the simulation was started with 
a nutrient saturated water column which initially gives 
rise to a uniform phytoplankton profile in the UML. 
Even though for the given parameter range the DCM 
is globally attracting, the center of biomass at first 
is trapped inside the UML as described above. This 
configuration can be sustained for a rather long du- 
ration. However, with ongoing nutrient depletion the 
phytoplankton density is slowly declining and as soon 
as the biomass in the UML is not suflftcient to shade 
light below it, the system undergoes a rapid transition 
to a DCM. Similar dynamics at two very diff'erent time 
scales are observed in Fig.[2f, where a DCM gradually 
moves upward until it reaches the bottom of the UML 
and then the biomass rapidly shifts into the upper layer. 
This transition from one solution type to the other can be 
very fast. It occurs on the biological time scale rg = yu^ ^ 
and takes approximately 10-50 days for the model in our 
parameter range. By contrast, the slow unstable dynam- 
ics before the transition is determined by the diff'usive 
time scale tb = (Zb - Z^ix)^ /(2Dd), which theoretically 
can last several years in deep waters. Obviously, in sys- 
tems with essential seasonal variability this transition 
will never be reached and the observed distribution will 
depend on the initial distribution of resources. 

The bottom panel in Fig. [2] shows oscillatory phyto- 
plankton profiles in a system with UML. If the oscil- 
latory state is deep enough then they are not aff'ected 
by the presence of a UML (compare Figs. [2^ and[2f). 
Interestingly however, the UML can support these os- 
cillations even if Do < Dmin, i.e. in a system where 
the population would die out without a UML (Fig. [2^1). 
The reason can be found in the fact that locally, within 
the well-mixed UML, diff'usivity is much larger than the 
persistence threshold, Du » Dmin- Now at the begin of 
the cycle, with practically no biomass the nutrient can 
freely diff'use towards the UML, where the population 
can outgrow sinking as soon as the nutrient concentra- 
tion has reached a critical level. Then, if the biomass 
in the UML is not limited by light, the depletion of the 
nutrient shifts the production layer downwards into the 
weakly mixed water, where diff'usivity is below the per- 



sistence threshold. There the population sinks further 
and finally declines because of light limitation, so that a 
new portion of the nutrient can reach the UML and the 
cycle starts anew. In this way the presence of a localized 
region with strong mixing at the top of the water column 
triggers an "oscillatory pump" for the biomass, with the 
eff'ect that the persistence threshold of the sinking pop- 
ulation disappears. This is remarkable because most of 
the time in the cycle the bulk of the biomass remains 
located in the weakly mixed lower layers. 

To investigate the system behavior in a large range 
of parameters we performed simulations for 900 pairs 
of {Nb, I in), (Nb, Dd) and (Nb, v). The results are pre- 
sented in the stability diagrams Figs.[3]\ -WP- As shown 
in Fig.[3]\, large values of I in in general lead to a DCM 
and large Nb to a UCM, while for intermediate re- 
source levels we observe a region with bistable behav- 
ior. The bistability range is reduced for smaller values 
of lin or Nb and disappears at a critical point (//„ ^ 
350 yumol photons m"^ s"^ and Nb ~ 25 mmol/m^ in 
Fig. [3K). Note that the lower border of the bistabil- 
ity range is quite well described by the analytical cri- 
terion Eq. (O, which is derived from the condition that 
the phytoplankton net production rate below the mixed 
layer is not positive (dash-dotted lines in Fig. [3]). In- 
terestingly, outside the region of bistability close to the 
critical point smooth transitions from deep to surface 
maxima are possible. In this case intermediate den- 
sity profiles with no clear separation between DCM and 
UCM appear and the phytoplankton biomass is located 
in both parts of the water column. To visualize the dis- 
tribution of biomass in this regime, parameter combina- 
tions, for which the median of the biomass distribution 
is located at Z^ix (i.e., half of the biomass is distributed 
within and half below the UML), are indicated as thick 
solid lines in Fig.O 

Fig. [3^ demonstrates the bistability range in the 
(Nb, Dd) parameter plane, where the bifurcation lines 
have an almost hyperbolic form. A UCM appears for 
large values of Do and Nb, whereas small values fa- 
vor a DCM. Very small values of Do result in oscilla- 
tory DCM (ODCM) solutions, in accord to the results 
of Huisman et al. (2006), and due to the presence of a 
UML this behavior can still be observed, and the popu- 
lation does not go extinct, for Do < Dmin- 

In the (Nb, v) parameter plane the bifurcation lines 
have almost vertical structure (Fig. [3]C). Most notably, 
the minimal value of Nb which can lead to upper 
biomass maxima practically is independent from v. 
Thus sinking has no strong influence on the stationary 
solutions in a well mixed layer. By contrast, oscillatory 
solutions arise only for a large enough sinking veloc- 
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ity. Again, the population can persist beyond the per- 
sistence threshold v > Vmax, where v^ax is the maximal 
sinking velocity which allows the survival of a popula- 
tion in a system without a UML. Vmax can be obtained 
from Eq. ([8]) for a given diffusivity. 

4. The two species model 

We now turn to the influence of a UML on the com- 
petition of two species which are difl'erently adapted to 
the conditions at the surface or in the deep water. For 
this, we extend the model to contain two phytoplankton 
species which difl'er in their respective half saturation 
constants Hi and H^ in Eq. ([3]). More specifically, we 
consider the competition between a most light-limited 
species (/-species) characterized by a low H^ and a 
large Hj value and so is forced to live close to the light 
source, and a most nutrient-limited species (A/^-species) 
with low Hi and large H^ that will usually do better in 
deep water. Note that the results of the previous section 
hold for the A/^- species. 

To avoid the influence of difl'erent growth, mortality 
and consumption rates, we study the simplest possible 
case and keep the other species' parameters identical in 
both species (see Table [T]). Thereby, in a well mixed 
uniform environment (e.g., a chemostat) these species 
would have parallel consumption vectors, so that the 
success of one species over the other depends only on 
the resource concentrations in the absence of biomass; 
with the consequence that in a chemostat such species 
cannot coexist and the outcome of their competition 
cannot be bistable (Tilman 1980, 1982). In a well mixed 
water column, where the light intensity reduces with 
depth, the competition of these species is more com- 
plicated, however the same results still hold (Huisman 
and Weissing 1995); whereas in an incompletely mixed 
water column, assuming only the light gradient, both 
species can coexist in a narrow parameter region (Huis- 
man et al. 1999). In the following we show that the 
outcome of competition is drastically altered in the pres- 
ence of two opposing resource gradients combined with 
a space dependent difl'usivity, resulting in new regions 
of competitive exclusion, bistability and coexistence. 

Figs. [3[) - [3^ present the two-species stability dia- 
grams in the difl'erent parameter planes, listing all possi- 
ble model outcomes. Basically, the overall shape of the 
transition lines between bistable and UCM/DCM states 
remain identical to those of the one-species model. 
The most notable difl'erence is that the bistability range 
(regime V) has become much wider and is extended 
toward smaller values of Nb (compare top and bottom 
panels in Fig.O. By contrast, the transition line to the 



UCM remains largely unchanged, since it is mainly de- 
termined by the minimal depth of the deep maximum 
at which it is not afl'ected by the UML. For our set of 
parameters, the deep maximum is always formed by the 
A/^- species (whereas the upper maximum can be formed 
by both species) and so this boundary does not change 
with addition of the /-species. Only in Fig. [3^ this 
boundary is essentially altered, reflecting the fact that 
sinking influences much stronger the species which oc- 
cupies the deep weakly mixed layers. 

More interestingly, the two species system exhibits 
a variety of patterns and new dynamical regimes (enu- 
merated by roman numerals I - VI in Figs. [3^ and[3f). 
A detailed representation of the biomass dynamics and 
phytoplankton profiles in each dynamic regime can be 
found in Fig.|4l 

These transitions and the species composition can be 
visualized in more detail by tracing the depth of the cen- 
ter of biomass as a function of Do- This is shown in 
Fig. \5\ for the single species and the two species case. 
This figure adds a third dimension to the vertical cross 
section through Figs. [3^ and [3^ at Nb = 20 mmol/m^. 
While the bulk of the biomass of the /-species mono- 
culture is always located within the UML (Fig. [5j\), a 
monoculture of the A/^- species for decreasing values of 
Dd undergoes four dynamical regimes, as described in 
the previous section, namely: UCM, UCM/DCM bista- 
bility, stationary DCM and ODCM (Fig. [5^). 

By contrast, the competition of two species leads to 
more intricate behavior, as is demonstrated in Fig. [SjC. 
At the high end of the Do range both species are located 
in the UML and can either coexist (region VI, Figs.|4l\ 
and |4^) or the A/^- species will competitively exclude 
the other species if light is the only limiting factor (not 
shown). With reduction of Do, the A^- species can form 
a DCM in the lower layers, yielding a large bistability 
range between a DCM of the A/^- species and a UCM of 
the /-species (region V, Figs.|4]C and|4t)). Thereby, the 
bistability of the phytoplankton profiles goes together 
with a bistability in the competition outcome. Follow- 
ing reduction of Do leads to competitive exclusion of 
the /-species and only stable DCMs of the A/^- species 
are found (region IV, Fig.|4^). 

For even smaller values of Dd the stationary DCM 
loses its stability yielding an oscillatory DCM (re- 
gion III, Fig. |4f ). Interestingly, in this regime, the /- 
species may obtain a time window during each cycle to 
establish a population in the UML (slightly visible as 
light blue stripes in Fig.|4f ), but the next rapid outburst 
of the A/^- species will again lead to the dominance of 
the A/^- species. With further reduction of Do, the period 
of the DCM oscillations increases until finally the N- 
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species cannot outcompete an established population of 
the /-species in the upper layer. However the oscillatory 
solution, once established, can persist for the same pa- 
rameter values. Thereby, we observe a second bistabil- 
ity regime, of either a stationary UCM of the /-species 
or the coexistence of both species due to oscillations (re- 
gion II, Figs.gt andllH). 

Finally reducing DjyXo very small values leads to a 
surprising result: the /-species always wins the com- 
petition and the steady UCM formed by the /-species 
becomes the single possible attractor in the system, in- 
dependent of the initial conditions (region I, Fig. O). 
Thus, a strongly light limited species is able to establish 
a steady population in the UML, outcompeting a less 
light limited species which, in the absence of the for- 



mer, would form an oscillatory maximum in the deep 
weakly mixed layers. Note that the /-species survives 
even iiDj) < Dmin, provided that the nutrient concen- 
tration in the upper layer is sufficiently high (Fig.|4K). 

Fig.[6]shows the outcome of competition between two 
species in the (A^g, /)/)) coordinate plane. In the system 
without a UML (Fig. [6j\) the A/^- species either wins or 
dominates (more than 75% of the biomass). Further- 
more, we observe no phytoplankton for diffusivity be- 
low Dmin- Comparison of Fig. [6]\ with Figs. [6^ and 
[6^! show that a UML plays a crucial role improving the 
competitive abilities of the /-species, which occupies 
a shallower depth because of its higher light require- 
ments. Shifts in the competition outcome can be caused 
by a change in the resource availability or by a bifur- 
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Figure 4: Typical phytoplankton profiles in the two species model in the presence of a UML. Shown are the time evolution of the total phytoplankton 
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cation in the vertical biomass distribution. The thick 
Hnes in Fig. [6^ and[6t! indicate those transitions in the 
species composition which are caused by a change of 
vertical biomass profiles. These lines coincide with the 
bifurcation lines in Fig. [3^. In the presence of a UML 
the /-species can win under the following conditions: 
(i) when the diff'usivity in the lower layers is reduced 



(Dd < O.lcm^/s) and the A/^-species alone would ex- 
hibit oscillating behavior; (ii) in the range of bistabil- 
ity, where either a UCM of /-species or a DCM with 
a large fraction of A/^- species appears; (iii) in the UML 
for sufficiently strong Do and low Nb ( 15 mmol n/m^ 
in Fig. [6^), when the A/^- species solely would still oc- 
cupy the UML, but becomes a weaker competitor under 
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stronger nutrient limitation. We also note that change 
of the vertical profiles causes non-monotonous shifts in 
the species composition. 

5. Graphical approach to analyze the competition 
outcome 

In a spatial system the outcome of competition is hard 
to deduce from first principles. First, the production 
layers of diff'erently adapted species typically will be 
located at diff'erent depths, which tends to reduce the 
strength of interspecific competition. Second, the com- 
petition becomes indirect, since the occurrence of phy- 
toplankton at a certain depth can completely reshape the 
resource distributions in other parts of the water column. 
To obtain an insight into the mechanisms of competi- 
tion in such an environment, we suggest a graphical ap- 
proach, which is based on a well known method devel- 
oped by Tilman for mixed systems (1980, 1982). 

If all resources and the biomass are uniformly dis- 
tributed, the state of the system may be represented as 
a point in a multidimensional resource space (Fig.[7]\). 
As the biomass grows and takes-up resources this sys- 
tem state point moves in the resource plane (along the 
consumption vector) until it hits the zero net growth iso- 
cline, characterized by a balance between growth and 
loss processes. From the condition m = jdi(N, I) in equa- 
tion ([3]) we can find such values N* and I* so that the 
specific growth rate equals the mortality 



These values determine the location of the zero net 
growth isoclines (Fig. [7]). The isoclines divide the re- 
source plane into areas of positive and negative local 
population net growth and mark possible resource com- 
binations in equilibrium. 

In a spatially-extended system, where the resource 
values change with depth z, the state of the system 
can be represented as a curve in the resource space, 
which is parametrically determined by the resource val- 
ues (N(z),I(z)). In this sense, the notion of a sys- 
tem state point in a homogeneous system naturally ex- 
tends to a system state curve (SSC) in a spatial ex- 
plicit system. A special case is given by a well mixed 
water column (Fig. |7]\): while the nutrients are uni- 
formly distributed, the light intensity decreases expo- 
nentially with depth and the SSC reduces to a line seg- 
ment {N = const, lout < I < lin}, where lout is the light 
intensity at the bottom of the water column (Huisman 
andWeissing 1995). 



To simplify the discussion, in the following we al- 
ways focus on the system state curve (SSC) in equi- 
librium. Fig. |7]\ shows a typical simulation outcome 
for the SSC of a single-species population in an incom- 
pletely mixed water column. Note that we use logarith- 
mic scaling to magnify the location of the SSC close to 
the zero net growth isoclines. Due to the diff'usive mix- 
ing the SSC does not settle at the zero net growth iso- 
clines, as would be the case for a uniform distribution of 
the resources. Instead the SSC extends into the area of 
positive ("favorable range"), as well as into the area of 
negative population net growth ("unfavorable range"). 
To give a crude insight into the density variation along 
a SSC we indicate the central range of the curve which 
contains 90% of the biomass as a thick line. 

For the analysis of a multi-species system, we find 
it convenient to calculate the SSC independently for 
each species in the absence of all other species. So 
each species attains its own SSC which would indicate 
its single- species resource configuration in equilibrium. 
Fig. [7^ shows an apparent example, illustrating the case 
of competitive exclusion. In this example the SSC of 
the A/^- species lies below the null growth isocline of I- 
species. Thus, the A/^- species reduces the resource levels 
at all vertical positions to values which are smaller than 
the resource requirements of the /-species and so do not 
permit positive net-growth of the invading /-species. By 
contrast, the /-species' SSC is located mostly above the 
null growth isocline of the A/^- species and so allows a 
positive net-growth of the A^-species. As a consequence, 
in this example the /-species will always be excluded as 
it has larger resource requirements. 

Fig. [Tf! demonstrates an example of the coexistence 
of two species. In this figure both curves intersect in 
such a way, that the SSC of one species is below that 
of its competitor in an essential part of its favorable 
range, and vice versa. In these conditions both species 
can coexist because both are superior competitors at 
diff'erent depths. In general, however, if the SSC of 
one species lies above the zero null growth isocline of 
its competitor, there is no unambiguous answer to the 
question "whether or not the latter species can invade?" 
The possibility of invasion will depend on the principal 
eigenvalue of a reaction-diff'usion equation, character- 
izing the growth of the biomass (Cantrell and Cosner 
2001; Ryabov and Blasius 2008), a problem which un- 
fortunately can only be solved in some simple cases. 

As is shown in Fig.jTt) the analysis of SSCs can also 
help to gain insight into the role of a UML. Here we 
use the same set of parameters as in Fig. [7^, but in the 
presence of a UML. In this case the main part of the /- 
species biomass is located close to the surface, while the 
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the axes. Parameter values: 
(C) Hj = 30 jumol photons m~^ ; 



30 mmol/m^ v = 0. (B, D) Hj = 98 //mol photons m" 
^ = 0.013 mmol nutrient m~^ for the /-species. 



s , //// = 0.02 mmol nutrient m ^ for the /-species; 



A/^- species still forms a DCM. Hence the presence of the 
UML leaves the SSC of the A/^-species unaffected while 
the SSC of the /-species becomes steeper and is shifted 
in direction of smaller resource values. This change in 



the form of the /-species' SSC results in new intersec- 
tions of both system state curves with the result to im- 
prove the competing abilities of the /-species. While the 
presence of the A/^- species still completely prevents an 
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invasion of the /-species, as in the case without a UML, 
now the SSC of the /-species Hes below the main part of 
the A/^-species' SSC. In this example, the A/^-species can 
exist in the presence of the /-species, however its frac- 
tion is small so that its contribution to the total resource 
distribution is negligible. As a result, in the presence 
of a UML the outcome of competition is bistable and 
depends on initial conditions. 

Our investigations show that the analysis of a species' 
SSC can further our understanding about resource com- 
petition in heterogeneous multi-species systems. What 
is more, one can obtain rigorous results for biologically 
similar species. A more detailed analysis will be pre- 
sented in a forthcoming article (Ryabov and Blasius, in 
preparation). 

6. Discussion 

In this article we investigated the influence of an up- 
per mixed layer on the distribution and competition of 
phytoplankton species in a water column, in which in- 
verse resource gradients (of light and a nutrient) can 
limit growth of the biomass. In this system the loca- 
tion of a production layer is not fixed, rather it depends 
on the initial and boundary conditions, on the species 
requirements, and on the stage of the process. Together 
with the presence of diff'erently mixed areas (i.e., the 
UML and lower layers), this leads to a plethora of phe- 
nomena, including bistability of phytoplankton profiles, 
changes in the competition outcome, and new critical 
conditions for the survival of a phytoplankton popula- 
tion. Thus our study not only proves the UML to be 
an important factor with the potential to shape the spa- 
tial distribution and species composition of phytoplank- 
ton, but also reveals insights of general ecological im- 
portance. 

While previous theoretical investigations have usu- 
ally focused on either a fully mixed or an incompletely 
mixed system, the presence of a UML requires a com- 
bination of these approaches. Including both factors, by 
dividing the water column into two separate compart- 
ments, Yoshiyama and Nakajima (2002, 2006) showed 
the existence of bistability in the spatial distribution of a 
phytoplankton monoculture. Our work confirms and ex- 
tends these findings on the following key points. First, 
assuming a gradual change of diff'usivity with depth, our 
modeling approach integrates the whole water column 
into a single framework. This allows us, for example, to 
investigate the influence of mixing in both the upper and 
the lower layers. Secondly, we analyze the competition 
of species which are difl'erently adapted to the availabil- 
ity of nutrients and light. As we show, this has drastic 



efl'ects because the species composition strongly corre- 
lates with the spatial patterning. And finally, our anal- 
ysis includes the case of a stratified lower layer, when 
a stable DCM cannot persist and oscillatory or chaotic 
solutions appear, which again have a strong influence on 
the species competition in the system. 

In comparison with lower layers three factors make a 
UML more favorable for phytoplankton species. Firstly, 
a UML has only one border with an unfavorable envi- 
ronment below the euphotic zone. In contrast, a deep 
production layer has two such boundaries and difl'usion 
of cells upward and downward from it leads to addi- 
tional losses due to either light or nutrient limitations. 
Secondly, the strong mixing in the upper layer allows 
for the survival of a sinking phytoplankton population. 
This additionally supports a population of phytoplank- 
ton even if the difl'usivity in the lower layers becomes 
very small. Higher difl'usivity also reduces the efl'ective 
mortality rate of sinking phytoplankton species (see Ap- 
pendix lA3]) . Thirdly, a UML promotes a nearly uniform 
distribution of nutrients, which makes the nutrient con- 
sumption more efficient and gives an additional com- 
petitive advantage for a species inhabiting a UML. All 
these factors decrease the total loss rate and, hence, re- 
duce the resource requirements. However, we note that 
a deep and turbid UML can also play a negative role, 
as it can lead to the extinction of species due to light 
limitation (Huisman et al. 1999b). 

Inclusion of a UML in the single species model al- 
lows for the existence of two alternative density config- 
urations with either a deep or an upper phytoplankton 
maximum. Note that in a small parameter range we ob- 
serve a third intermediate profile with an essential part 
of the biomass above and below the thermocline. More- 
over, there is a range of parameters for which the system 
is bistable and the appearence of either deep or surface 
biomass maximum is determined by the initial distribu- 
tion of nutrients. Both the consumption of nutrient and 
the self-shading of light are necessary conditions for this 
behavior. Since the biomass obstructs the upward nutri- 
ent flow, it makes the upper layer unfavorable, provided 
that a deep maximum of phytoplankton has established. 
The shading of light is an opposite mechanism, which 
prevents population growth in the lower layers. The 
third requirement for bistability is the presence of strong 
mixing in the upper layer. This decouples the locations 
of the production layer from that of the bulk biomass 
and prevents a drift of the population toward a DCM, 
as it would occur in a system without a UML. Similar 
bistable behaviour has also been suggested to occur in 
field data. For example, in a recent study of Antarctic 
waters (Holm-Hansen and Hewes 2004) it was found 
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that adjacent spots could randomly exhibit surface max- 
ima or DCM, even though temperature, salinity and wa- 
ter density profiles of the water columns were practi- 
cally identical. 

The occurrence of bistability may have important 
ecological consequences and usually goes together with 
catastrophic shifts. Imagine that we slowly change a 
parameter, i.e., the incident light intensity, so that its 
value passes through the bistability range. Then there 
will be a point, where the current state suddenly loses 
its stability, thereby inducing a rapid transition from one 
phytoplankton profile to another. This is also impor- 
tant for a border between oligotrophic and eutrophic wa- 
ters, where this transition might be induced by seasonal 
changes and leads to a shift of the interface between a 
deep and a surface biomass. maximum. In this scenario 
bistability will lead to a lag in the transition from one 
state to another and back. 

Outside the bistability range the system possesses 
only one attractor. However, the transition from an un- 
stable to a stable solution may take a long time (espe- 
cially close to the bistability range) and includes two 
stages. During the first stage, the nutrient concentration 
slowly evolves to a level which allows for the forma- 
tion of a stable biomass profile. This stage occurs on 
the slow diff'usive time scale and may typically last 5- 
10 years in deep waters. Once the nutrient values have 
crossed a critical value, a rapid transition to the sta- 
ble profile is triggered. This second stage develops on 
the biological time scale and lasts approximately 10-50 
days. 

As a consequence, in this, and in the bistable, regime 
the system is very sensitive to disturbance. If, due to 
some perturbation, a large quantity of nutrients is trans- 
ported into the UML, the system can rapidly switch 
from a DCM to a UCM and remain in this state for a 
long time. Similar, a prolonged lack of light can induce 
a UCM, when the reduced phytoplankton density can 
not prevent diflfusion of nutrients to the surface. While 
outside of the bistability range the reverse transition in 
principle is possible without further external perturba- 
tion, such a flip back to the original state will usually 
take a long time. Thus, the eflfective size of the bista- 
bility regime, i.e. the parameter regime in which two 
dynamic states are preserved for long times and dis- 
turbance can induce long-lasting shifts in the system, 
is much larger than the true bistability range and com- 
prises a large part of the model's parameter space. Since 
the relaxation time is long, in a seasonal environment 
such a system might exhibit only an unstable state (most 
probable an upper chlorophyll maximum). 

The eflfects of a UML are even more pronounced in 



a system that includes two competing species which are 
diff'erently adapted to light and nutrient limitation. To 
analyze this situation we presented a graphical approach 
which helps to understand the mechanisms of competi- 
tion in spatially extended systems. Moreover, we found 
that in the range of parameters where the two species 
can independently form an upper or a deep biomass 
maximum, the two species model demonstrates bista- 
bility both in the spatial distribution and in the compe- 
tition outcome. Compared to the single-species model, 
the bistability range is considerably enlarged. Note that 
the bistability in the competition outcome is induced 
by the UML, whereas in a homogeneously mixed water 
column these species either coexist or only one species 
wins. 

Another remarkable finding is the survival of a sink- 
ing phytoplankton population even if the diff'usivity 
in the deep layers cannot prevent population washout 
(Shigesada and Okubo 1981; Speirs and Gurney 2001; 
Straube and Pikovsky 2007). While the retention of 
phytoplankton in the upper layer can be provided by 
many factors such as the regulation of algal buoyancy or 
the formation of eddies, etc (Raymont, 1980), our model 
shows another dynamic mechanism: in the absence of 
phytoplankton in the deep layers, the nutrient can dif- 
fuse upward into the UML, where the population can 
start to grow as the nutrient concentration reaches a suf- 
ficient level. However if the light limitation is not strong 
enough the biomass shifts into the deeper layers where it 
cannot outgrow the sinking and ultimately declines be- 
ing limited by light. Then the nutrient again can diff'use 
upwards and the cycle repeats. Thus we identify a novel 
mechanism of how a population can overcome the drift 
paradox, as the persistence threshold disappears in the 
presence of a UML. 

The mechanisms ensuring persistence ultimately are 
related to the fact that locally, within the well-mixed up- 
per layer, diflfusivity is much larger than the persistence 
threshold. Thus obviously, a stationary distribution of 
a strongly light limited species within the UML will 
not be affected by low values of diff'usivity in deeper 
layers. However more astonishingly, that the persis- 
tence threshold is absent for the oscillatory solutions 
for which the bulk of the biomass remains located in 
the weakly mixed lower layers for most of the time 
in the cycle. Again, this eff'ect is related to the mo- 
bility of the production layer, which for an oscillatory 
DCM undergoes cycles in depth. These vertical cycles 
of the production layer can act as an "oscillatory pump", 
so that during some, possibly very small, phase of the 
cycle the production layer is located within the UML. 
These short cyclic visits of the UML are suflficient to 
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refuel biomass and to prevent extinction of the popula- 
tion. Thereby, however, the biomass can be temporarily 
reduced to very small values so that different model out- 
comes might be expected by including an Allee effect. 

Finally, we have identified an interesting competitive 
exclusion effect, where a species can be outcompeted 
in dependence of its dynamical state. This occurs for 
low values of deep diffusivity, when the DCM of the A^- 
species becomes oscillatory (Huisman et al. 2006) and a 
further reduction of diffusivity favors to the domination 
of the strongly light limited /-species, which still occu- 
pies the UML. At the first glance this result is counterin- 
tuitive, because usually a decrease of nutrient transport 
results in an increased depth of the biomass maximum 
and here we suddenly obtain a UCM solution instead 
of a DCM. The replacement of the oscillatory solutions 
occurs because the /-species forms a non-oscillatory 
UCM and shades light, inhibiting the rapid outburst of 
the biomass of the A/^-species in the lower layers. As a 
consequence the (non-oscillatory) UCM configuration 
in the UML is able to replace the oscillatory deep max- 
imum. As some climate models predict higher water 
stratification with an increase of temperature (Bopp et 
al. 2001; Sarmiento et al. 2004), our findings may help 
to predict and understand future changes in phytoplank- 
ton patterns ongoing with global climate change. 

Possible interesting further extensions of our model 
would be the inclusion of zooplankton or higher trophic 
levels, the inclusion of an Allee effect or the extension 
of our study into a three-dimensional system. 
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After integration of these equations with respect to the 
boundary conditions, Eqs. Q and d?]). 
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we obtain the total biomass 



P(z)dz = 



Dr 



am(l - e) dz 



z^Zb 



(12) 



This formula relates the phytoplankton biomass and the 
nutrient flux and can be interpreted as a conservation 
law in the system. 

To describe the nutrient flow we first consider the re- 
gion, which is free of the phytoplankton biomass. Using 
Eq. ([TT1) we obtain 



J-^yz) — - — = const . 
oz 



(13) 



This equation corresponds to the stationary diffusion 
flow and goes together with a linear gradient of nutri- 
ents, below the phytoplankton maximum. If the total 
phytoplankton biomass is located within the UML, the 
nutrient flow can be estimated as 



D{z) 



dN(z) 



Dr 



Nb-N* 



(14) 



dz Zb — Zjnix 

where N* from Eq. (O approximates the concentration 
of the nutrient within the UML. Substituting this into 
Eq. ([T2I) , we obtain the total phytoplankton biomass in 
the UML (see Fig.©. 



W = 



Dr 



am(l - s) Zb - Z^ 



(15) 



A. Analytical derivations 

A.L Total equilibrium phytoplankton biomass in the 
single-species model 

At equilibrium the left hand sides of Eqs. ([B and Q 
equal zero 



dP d 

li{NJ)P-mP-v— + — 
oz oz 



-a/i(N, I)P -h samP -h — 
oz 



dP 
D(z) — 
dz 

D(z)— 
oz 



= 0(10) 



= 0(11) 



A. 2. Border of the stability of a UCM 

A sufficient condition for stability of a UCM is the 
light limitation of growth below the UML, i.e., the light 
intensity below the UML should be smaller than the crit- 
ical light intensity /*, see (|9]). Using Eq. (|4]), we obtain 



,exp 
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(16) 



If the total phytoplankton biomass is located in the 
UML, equations ([T5]) and ([T6l) give the following cri- 
terion for stability of the upper maximum 
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Figure 8: Total phytoplankton biomass W as a function of the bot- 
tom nutrient concentration A^^ for different values of the deep turbu- 
lent diffusivity Do- Comparison of the results from (T5\ (solid lines) 
with numerical simulations (symbols). The analytic estimation is in 
excellent agreement with the numerical results with exception of the 
region of high diffusivity Do and nutrient concentration A^^, in which 
the growth of phytoplankton biomass is limited by light and does not 
depend on the nutrient concentration. Thus, the concentration of nu- 
trients in the UML will be larger than N* and the estimation ( flU for 
the nutrient gradient fails. 

This line is shown in Figs.[3]\ -[3f as the lower bound- 
ary of the bistabihty range and is in a good agreement 
with the numerical simulations. 

A. 3. Losses in the UML 

The upper mixed layer is more favorable for sinking 
phytoplankton species. To show this, consider a water 
column where the diffusivity does not depend on depth. 
Its is easy to show that washout from the water column 
can be interpreted as an additional mortality term. Sub- 
stituting P = Pexp (vz/2D) into Eq. ([T]), we obtain 

dP ~ / vM ~ d^P 

Note that dtP has the same sign as dtP, thereby both 
functions grow and decline simultaneously. Introduce 
the new mortality 

m' =m + — , (18) 

and substitute it in the expression for the limiting re- 
source values ©. Since m' > m, the new limiting val- 
ues /'* and N'* in the presence of sedimentation should 
be larger, see (|9]). This results in a shift of the zero 
net growth isoclines towards higher values of resources. 



However, in a well mixed layer the term IAD van- 
ishes, leading to lower resource requirements. Thus a 
UML creates more favorable conditions for the sinking 
phytoplankton biomass. 
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Table 1: Parameters values and their meaning 



Symbol 


Interpretation 


Units 


Value 




Independent variables 








t 


Time 


h 


- 




z 


Depth 


m 


- 




Dependent variables 








P(zJ) 


Population density 


cells m~^ 






KzJ) 


Light intensity 


//mol photons m"^ s"^ 






N(z, t) 


Nutrient concentration 


mmol nutrient m"^ 






Parameters 








/in 


Incident light intensity 


jimol photons m~^ s~^ 


600(100- 


600) 




Background turbidity 


m-i 


0.045 




k 


Absorption coefficient of phytoplankton 


m^ cell"^ 


6x10-10 




Zb 


Depth of the water column 


m 


300 






Depth of the upper mixed layer 


m 


50 




W 


Characteristic width of the thermocline 


m 


1 




Dd 


Turbulent ditfusivity in the deep layers 


? -1 
cm s 


0.3 (0.04 - 


1) 


Du 


Turbulent diffusivity in the UML 


2 -1 

cm^ s ^ 


50 




/^max 


Maximum specific growth rate 


h-i 


0.04 




Hi 


Half saturation constant of light 


jumol photons m~^ s~^ 


20; 98 






limited growth for N- and /-species 










Half saturation constant of nutrient 


mmol nutrient m~^ 


0.0425; 0.015 




limited growth for A^- and /-species 








m 


Specific loss rate 


h-i 


0.01 




a 


Nutrient content of phytoplankton 


mmol nutrient cell"^ 


1 xlO-9 




s 


Nutrient recycling coefficient 


dimensionless 


0.5 




V 


Sinking velocity 


mh-i 


0.042 




Nb 


Nutrient concentration at Zb 


mmol nutrient m"^ 


5-100 
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Table 2: Acronyms 


Symbol 


Interpretation 


DCM 


Deep chlorophyll maximum 


ODCM 


Oscillatory or chaotic 




deep chlorophyll maximum 


UCM 


Upper chlorophyll maximum 


UML 


Upper mixed layer 


SSC 


System state curve 
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